      Subroutine BORROW_FREEFILE(FREENUM)

      include "param.fi"
      include "common.fi"

      INTEGER I, FREENUM
      I = 1
      Do While (FREESIGN(I) .EQ. 1)
        I = I + 1
      Enddo

      FREENUM = I

      FREESIGN(I) = 1

      End Subroutine

C==================================================================
      Subroutine RETURN_FREEFILE(FREENUM)

      include "param.fi"
      include "common.fi"

      Integer FREENUM

      FREESIGN(FREENUM) = 0

      End Subroutine
C==================================================================
      Function HASHIND(IND)
      ! Periodic condition
      include "param.fi"
      include "common.fi"
      INTEGER IND
      IF (IND.LE.0) Then
        HASHIND = IND + GRIDNUM
      ELSEIF (IND.GT.GRIDNUM) Then
        HASHIND = IND - GRIDNUM
      ELSE
        HASHIND = IND
      Endif
      End Function
C==================================================================
      Function HASHIND_Ghost(IND)
      ! Ghost particle condition
      include "param.fi"
      include "common.fi"
      INTEGER IND
      IF (IND.LE.0) Then
        HASHIND_Ghost = IND - 1
      ELSEIF (IND.GT.GRIDNUM) Then
        HASHIND_Ghost = - (GRIDNUM + GRIDNUM + 1 - IND)
      ELSE
        HASHIND_Ghost = IND
      Endif
      End Function
C==================================================================
      Subroutine HASHINDEX31(I1,I2,I3, I123)
      include "param.fi"
      include "common.fi"
      Integer I1, I2, I3, I123
      I123 = (I3-1)*GRIDNUM*GRIDNUM + (I2-1)*GRIDNUM + I1
      End Subroutine
C==================================================================
      Subroutine HASHINDEX13(I123, I1,I2,I3)
      include "param.fi"
      include "common.fi"
      Integer I1, I2, I3, I123, IJK
      I1 = MOD(I123, GRIDNUM)
      If (I1.EQ.0) I1 = GRIDNUM
      IJK = (I123 - I1) / GRIDNUM
      I2 = MOD(IJK, GRIDNUM) + 1
      IJK = (IJK - I2 + 1) / GRIDNUM
      I3 = IJK + 1
      End Subroutine
C==================================================================
      Subroutine NP2DINDEX31(I1,I2,I3, I123)
      include "param.fi"
      include "common.fi"
      Integer I1, I2, I3, I123
      I123 = (I3-1)*NP + I1
      End Subroutine
C==================================================================
      Subroutine NP2DINDEX13(I123, I1,I2,I3)
      include "param.fi"
      include "common.fi"
      Integer I1, I2, I3, I123, IJK
      I1 = MOD(I123, NP)
      If (I1.EQ.0) I1 = NP
      IJK = (I123 - I1) / NP
      I3 = IJK + 1
      I2 = 0
      End Subroutine
C==================================================================
      Function Distance2P(XA, YA, ZA, XB, YB, ZB)
      include "param.fi"
      include "common.fi"
      REAL XA, YA, ZA, XB, YB, ZB, XX, YY, ZZ
      XX = ABS(XA - XB)
      If (XX.GT.HALFCUBESIZE) XX = XX - CUBESIZE
      If (XX.LT.-HALFCUBESIZE) XX = XX + CUBESIZE
      YY = ABS(YA - YB)
      If (YY.GT.HALFCUBESIZE) YY = YY - CUBESIZE
      If (YY.LT.-HALFCUBESIZE) YY = YY + CUBESIZE
      ZZ = ABS(ZA - ZB)
      If (ZZ.GT.HALFCUBESIZE) ZZ = ZZ - CUBESIZE
      If (ZZ.LT.-HALFCUBESIZE) ZZ = ZZ + CUBESIZE
      Distance2P = SQRT(XX*XX+YY*YY+ZZ*ZZ)
      End Function
C==================================================================
      Subroutine Distance2P_VEC(XA, YA, ZA, XB, YB, ZB, XX, YY, ZZ)
      include "param.fi"
      include "common.fi"
      REAL XA, YA, ZA, XB, YB, ZB, XX, YY, ZZ
      XX = XA - XB
      If (XX.GT.HALFCUBESIZE) XX = XX - CUBESIZE
      If (XX.LT.-HALFCUBESIZE) XX = XX + CUBESIZE
      YY = YA - YB
      If (YY.GT.HALFCUBESIZE) YY = YY - CUBESIZE
      If (YY.LT.-HALFCUBESIZE) YY = YY + CUBESIZE
      ZZ = ZA - ZB
      If (ZZ.GT.HALFCUBESIZE) ZZ = ZZ - CUBESIZE
      If (ZZ.LT.-HALFCUBESIZE) ZZ = ZZ + CUBESIZE
      End Subroutine

C==================================================================
      Subroutine Distance2P_VEC_Ghost(XA, YA, ZA, XB, YB, ZB,
     &                                GI, GJ, GK, XX, YY, ZZ)
      include "param.fi"
      include "common.fi"
      REAL XA, YA, ZA, XB, YB, ZB, GI, GJ, GK, XX, YY, ZZ
      REAL XXA, YYA, ZZA
      XXA = XA
      YYA = YA
      ZZA = ZA
      ! Ghost operation for location
      If (GI.LT.0.0) Then
        If (XXA.LT.HALFCUBESIZE) XXA = - XXA
        If (XXA.GE.HALFCUBESIZE) XXA = CUBESIZE + CUBESIZE - XXA
      Endif
      If (GJ.LT.0.0) Then
        If (YYA.LT.HALFCUBESIZE) YYA = - YYA
        If (YYA.GE.HALFCUBESIZE) YYA = CUBESIZE + CUBESIZE - YYA
      Endif
      If (GK.LT.0.0) Then
        If (ZZA.LT.HALFCUBESIZE) ZZA = - ZZA
        If (ZZA.GE.HALFCUBESIZE) ZZA = CUBESIZE + CUBESIZE - ZZA
      Endif
      XX = XXA - XB
      YY = YYA - YB
      ZZ = ZZA - ZB
      End Subroutine
C==================================================================
      Function Dual2to1(I, J)
      include "param.fi"
      include "common.fi"
      Integer I, J
      Dual2to1 = (J - 1) * GridDual + I
      End Function
C==================================================================
      Function PeriodMid(XX1, XX2)
      include "param.fi"
      include "common.fi"
      Integer I, J
      REAL XX1, XX2, XXMIN, XXMAX

      If (XX1 .LT. XX2) Then
        XXMIN = XX1
        XXMAX = XX2
      Else
        XXMIN = XX2
        XXMAX = XX1
      Endif

      If ((XXMAX-XXMIN) .LT. (XXMIN-XXMAX+CUBESIZE)) Then
        PeriodMid = 0.5 * (XXMIN + XXMAX)
      Else
        PeriodMid = 0.5 * (XXMIN + XXMAX + CUBESIZE)
        If (PeriodMid.GT.CUBESIZE) PeriodMid = PeriodMid - CUBESIZE
      Endif

      End Function
C==================================================================
      Function InBox(XX)
      include "param.fi"
      include "common.fi"
      Real XX

      InBox = XX
      If (InBox.LT.0.0) InBox = InBox + CUBESIZE
      If (InBox.GT.CUBESIZE) InBox = InBox - CUBESIZE

      End Function
